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Abstract. We discuss a method for constructing generalized Wannier functions that 
^ I are maximally localized at the minima of a one-dimensional periodic potential with a 

^ ■ double-well per unit cell. By following the approach of (Marzari M and Vanderbilt D 

1997 Phys. Rev. B 56, 12847), we consider a set of band-mixing Wannier functions 



with minimal spread, and design a specific two-step gauge transformation of the Bloch 
functions for a composite two band system. This method is suited to efhciently 
computing the tight-binding coefficients needed for mapping the continuous system 
■ to a discrete lattice model. Their behaviour is analyzed here as a function of the 

^ ■ symmetry properties of the double- well (including the possibility of parity-breaking), 

I in a range of feasible experimental parameters. 

(N 

> ; PACS numbers: 67.85.Hi, 03.75.Lm, 03.65.Vf 
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1. Introduction 

Ultracold atoms in optical lattices are attracting an increasing interest as quantum 
simulator for condensed matter systems [H Optical lattices are realized by means 
of one or more continuous sinusoidal potentials created by exploiting the electric dipole 
ij interactions between the atoms and the laser beams pQ. Depending on the beam 

^ I geometry one can realize one-, two-, or three-dimensional periodic lattices, with one 

or more wells per unit cell [3]; quasiperiodic structures are also possible (see e.g. 

[21 a El IS m). 

In the tight binding regime, when the lattice intensity is sufficiently high so that the 
atoms are deeply localized in the lowest vibrational states of the potential wells (each 
well being associated to a site of a discrete lattice), it is convenient to map the system 
hamiltonian onto a discrete lattice model (or tight binding model), representing the usual 
tool for theoretical calculations. The paradigms are the Hubbard model for fermions 
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[8], and the Bose-Hubbard model for bosons [9J. These models are characterized by 
tunnelling coefficients related to the hopping between neighbouring sites, and interaction 
strengths which characterize the onsite interaction among the atoms, whose actual values 
depend on the parameters of the underlying continuous model. 

As the mapping between the continuous and discrete versions of the system 
hamiltonian is achieved by means of an expansion over a basis of localized functions at 
each potential well PQ dUl [H] , a precise knowledge of these basis functions is therefore 
important to connect the actual experimental parameters with the coefficients of the 
discrete model employed in the theoretical calculations. 

For a simple sinusoidal potential, the natural basis is provided by the exponentially 
decaying Wannier functions discussed by Kohn [121 [13] • Notably, in this case the 
expression for the tunnelling coefficient turns out to depend just on the Bloch spectrum 
[T4j . being therefore independent on the basis choice (this does not hold for the 
interaction coupling). Analytic expressions for both coefficients can be obtained by 
means of different approximations fU] [T5]. 

In the case of two wells per unit cell, the Kohn- Wannier recipe is not sufficient. 
For example, for a symmetric double- well the Kohn- Wannier functions display the same 
symmetry as the local potential structure [131 [IB] and thus they occupy both wells and 
cannot be associated to a single lattice site. A common approach used in the literature 
is that of the so-called atomic orbitals [191 HZHH], that has been recently employed e.g. 
for the case of a symmetric double-well unit cell of two-dimensional graphene-like optical 
lattices [20]. This method is based on a specific ansatz, according to which tight-binding 
Wannier functions are constructed from linear combinations of wave functions deeply 
localized in the two potential wells of the unit cell. 

A more general approach is the one proposed by Marzari and Vanderbilt [21], 
where maximally localized Wannier functions (MLWFs) are obtained by minimizing the 
spread of a set of Wannier functions by means of a suitable gauge transformation of 
the Bloch eigenfunctions. This method coincides with the Kohn method for the single 
band MLWFs of a one-dimensional potential, but it can be extended to more complex 
situations when generalized MLWFs for composite bands are needed. The method is 
implemented by means of a software package, and is largely employed for computing 
MLWFs of real condensed matter systems [22] . 

In this article, we consider the case of a one- dimensional periodic potential with a 
double- well per unit cell [231 1211 ISSl [261 123 [2H1 [29] , discussing an alternative method for 
constructing the low-lying generalized MLWFs based on the minimal spread requirement 
of Marzari et al. [21], specifically suited for double- well potentials. In particular, we 
consider a composite band formed by the two lowest Bloch bands and, differently from 
j21j . we design a two-step gauge transformation specific for a composite two band system, 
that can be solved by integrating a set of ordinary differential equations, with suitable 
boundary conditions. This allows to efficiently compute the tunnelling coefficients and 
the other tight binding coefficients in terms of the parameter of the continuous potential. 
We also remark that, though the approach of Marzari et al. [21j was proposed for 
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degenerate bands, we find that the method works properly in a wider regime, provided 
that the first band gap does not exceed the distance between the second and third band. 

The paper is organized as follows. In Section |2] we review the mapping of a many- 
body hamiltonian onto a tight binding model, discussing in particular the case of a 
sinusoidal periodic potential with a double well in the unit cell. Then, in Section [3] we 
describe the gauge transformation for constructing the generalized MLWFs, giving some 
examples and comparing them with those of the single band approach. Then, in Section 
m we discuss the regime of validity of the composite band approach, and compare the 
predictions of the full and nearest-neighbour versions of the model by discussing the 
behaviour of the tunnelling coefficients and the other tight-binding parameters. Final 
considerations are drawn in the Conclusions. Technical points regarding the mapping on 
momentum space and the numerical implementation are discussed in the Appendices. 

2. Tight binding models for double well optical lattices 

Let us start by reviewing how tight-binding models are defined from a continuous 
potential. As a specific example, here we consider a one- dimensional many-body 
hamiltonian for ultracold bosons [10] 

u = jdx ^^Ho^ + 1 /^^ V'^V'^^^ = no + n,at (1) 

with iIj{x) being the bosonic field operator, Hq = —{h"^ /2m)'V^ + V{x) the single particle 
hamiltonian, and V{x) a double-well periodic potential of the form 

Vix) = Vi sin^ {kBX + 0o) + V2 sin^ {2kBX + 9o + 20o) (2) 

with ks = T^/d and Vi strictly non vanishing (Vi > 0) in order to fix the overall period 
to d, V{x + d) = V{x). This is a typical potential used in the experiments with ultracold 
atoms [T], by which one can construct an optical lattice with one or two wells in the unit 
cell. In particular, the condition for having two minima in each period for any value of 
^0 is V2 > 0.5Vi. In the following, the potential amplitudes Vi will be expressed in units 
of En = h'^k'^/2m, the so-called recoil energy for an atom absorbing a photon of the 
first lattice. 

The phases Oq and 0o are arbitrary; (pQ represents a rigid shift of the whole potential, 
while 6*0 can be varied along with the ratio of the amplitudes V2/V1 to change the 
landscape of the potential. For convenience, the unit cell is defined as having two 
maxima at the cell borders, with both minima inside the cell (see figure [1]). In addition, 
the angle 0o is tuned in order to have a unit cell centered in a; = 0, x G [—d/2, d/2], with 
the absolute minimum in the left well 0. Depending on the value of ^0, it is possible 
to realize three different configurations, as shown in figure [TJ (a) A unit cell with two 
different minima and with degenerate maxima, for = nvr (n G Z); in this case the 
whole periodic potential has two (classes of) parity centers, corresponding to the two 

I The presence of 0o is also useful for testing the robustness of the numerical method for obtaining the 
MLWFs, that indeed should not depend on an overall translation of the potential. 
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Figure 1. The three possible configurations for the unit cell of the potential in ([2]) 
(here V2 = 2Vi): (a) two different minima, with the overall potential having two centers 
of parity, for 9o — (0o — 7''/4); (b) an asymmetric double- well with parity that is 
broken globally - in this example 6*0 = 7r/4 (0o ~ 7r/8); (c) a symmetric double- well, 
Oq = 7r/2 ((^0 = 0). The black dots in (a), (c) represent the parity centers of the whole 
periodic potential. 

minima, with all the maxima being degenerate, (b) An asymmetric double-well with 
parity that is globally broken, for any ^0 ^ (0, 7r/2) -|-n7r/2. (c) A symmetric double- well 
in the unit cell, for 6q = 7^/2 + nvr; in this case the potential has again two centers of 
parity, now at the two maxima, with all the minima being degenerate. 

As anticipated in the Introduction, when the potential wells are deep enough, it 
may be convenient to map the hamiltonian ([1]) onto a tight binding model on the discrete 
lattice corresponding to the potential minima, by expanding the field operator ip{x) on 
a basis {fnj{x)} of functions localized around each minimum 

= ^dnjfnj{x) (3) 

nj 

where alj {o^nj) represent the creation (destruction) operator of a single particle at site 
j, and satisfy the usual commutation rules [Snj, = Sjj'6nn' (following from those 
for the field ■?/'). 

In the presence of a single well per unit cell, it is known that a basis of localized 
functions is provided by the exponentially decaying Wannier functions Wnj{x) discussed 
by Kohn [131 E] • In general, this is not the case when there are two wells per unit cell. 
For example, for a symmetric double- well the Kohn- Wannier functions display the same 
symmetry of the local potential structure [131 HB] and thus they cannot be associated 
to a single lattice site as as they occupy both wells in the unit cell. In the next section 
we will show that when the two lowest Bloch bands are sufficiently close to each other 
with respect to the third band (as it will be clear from the discussion in Sect. 13.21 and 
S]), we can construct a set of generalized Wannier functions Wnj{x) that are maximally 
localized at each minima, by following the approach of Marzari and Vanderbilt [21] for a 
composite band. This corresponds to the generalization of the single band approximation 
(in case of a single well lattice) to the double well case, as we need at least two localized 
functions in each lattice cell to map the system on the discrete lattice. Then, in Section 
m we will explore the range of validity of this composite band approaches, highlighting 
the different implications on the structure of different tight binding models. 
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Figure 2. A sketch of the double- well structure and of the tunnelling coefficients 
between sites A and B. 

In the following we will restrict the analysis to the two lowest energy bands, in 
analogy with the single band approximation for the Bose-Hubbard model [11]. Then, 
within this approximation, the single particle hamiltonian can be written as 

^0- Xl^jW^'^-^J^'l^ol/iv) (4) 

uu'=A,B jj' 

where j is the unit cell index whereas u = A,B substitutes the band index n = 1,2 
being an internal index labelling the left and right sub- wells respectively (see figure |2]). 
Here the expansion coefficients correspond to the onsite energies E,^ = (fj^Holfj^), and 
to the tunnelling amplitudes between different (sub) wells = —{fjy\Ho\fji^i). In 
general, it is customary to further approximate the above expression by neglecting the 
coupling beyond nearest neighbours both for the single- well [10] and double- well lattices 
[25| l28t [29] . This is a reasonable assumption for a single well lattice in the tight binding 
regime [HI |31], but may not be fully justified in the range of the typical experimental 
parameters for a double well, as we will discuss later on. For this reason, here we keep 
all the terms corresponding to nearest-neighbouring cells, characterized by the following 
tunnelling coefficients 

= - (/^.J^ol/a+D.) (5) 

Tab = -(/mI^oI/,.) (6) 

Jab± = - (/jA 1^0 1 /(j±i)fl) (7) 
as shown in figure |2j The corresponding hamiltonian is 

u=A,B j u=A,B j 
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(TABa]^ajg + JAB+«]^o,(j+i)s + -^AB_«j^a(j-i)s + h.c. 



When the tunnelling and Jab+ are negligible the above expression can be further 
simplified by retaining just the coupling between nearest neighbouring wells 

'^0- Y Y ~ Y (^^e^L^is + JAB.a]^a(^j_i)^ + h.c^ . (9) 

u=A,B j j 

This is the analogous of the nearest-neighbours approximation for the single-well case 
and it is commonly used in the literature [281 ES]. Hereinafter, we will refer to the 
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approximate hamiltonians in (|8]) and ([9]) as the (single particle) full and nearest- 
neighbour tight-binding models, respectively. 

Finally, the general form of the interaction term is 

'Hint = — 1^1 ^Li^2%3i^3%4'^4 / /j2!^2 /j3 1^3/34 1^4 • (-l-O) 

Here we consider just the onsite contributions in the two sub-wells A and B 

n.nt - J2 " ^^^^ 

u=A,B ju 

with Uy = g j dx\fj^{x)\^ . This corresponds to the usual approximation used for the 
single well case (see e.g. [lO]), and is justified when the basis functions {/nj(a^)} are 
strongly localized in each sub-well. In general, once these functions are known, the 
next-to-leading interaction terms can be straightforwardly computed by evaluating the 
corresponding superposition integral of four /^/s. 

2.1. Single particle spectrum 

Let us now consider the single particle spectrum of the hamiltonian (|8]). By defining 




27r 



the full tight binding hamiltonian (|8]) can be rewritten as 

no = y2 [ dk K,,{k)bibku' (13) 



with 



and 



e^{k) = E^-2J^cos{kd) (15) 

Z{k) ^ -{Tab + Jab^c-'"' + JAB.e'""), (16) 

where the operators b^u satisfy the canonical commutation relations, [buk, b^'k'] = ^{^ " 
k'^byyi. Then, by diagonalizing the matrix h{k), and defining e±(A;) = (e^(/c) ±eB(A;))/2, 
we get (see also [2(1]) 

el{k) = e^k) ± JeUk) + \Z{kW (17) 



that represents the spectrum of the full tight binding model in (|8]). In addition, with 
Jy = = Jab+, the same expression gives also the spectrum for the nearest-neighbour 
approximation in (Q. 

Instead, in the single band case we simply have [1] 

e:\k) = E:,' - 2J:!> cos{kd) (18) 
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with 



n 



isb 




dk en{k); 



J: 



n 



sb 




dk en{k)e' 



ika 



(19) 



en{k) being the exact Bloch spectrum; notably these expressions do not depend on the 
choice of the Wannier basis. 

3. Generalized Wannier functions 

In this section we discuss the method for constructing the MLWFs for the double well 
case. In order to fix the notations, let us first recall some basic properties of periodic 
systems [19l[30]. Owing to the Bloch's theorem, the eigenf unctions of the single particle 
hamiltonian Hq can be written as ipnk{x) = e^^^Unki^)-, the Unki^Y^ having the same 
periodicity of the potential and satisfying the following normalization in the unit cell, 
{umk\unk) = {d/2n)6mn- The Wannier functions for a single band are defined as 



with B indicating the first Brillouin zone, k G [— A;^,^^], and Rj = jd, whereas 
generalized Wannier functions for composite bands have the same formal expression 
but are built up from a linear combination of Bloch eigenstates, namely 



with UW = 1 and Unmik + 2kB) = Unmik). The Wannier functions satisfy the 
ortho-normality relation (w„j|w„'j') = {wnj\wn'ji) = Snn'^jj'- We also remark that the 
generalized Bloch functions ipnk are not eigenstates oi Hq; however, their ortho-normality 
relations are preserved, owing to the unitarity of the transformation matrices Unm{k)- 

Different choices of the matrices Unm{k) lead to different results and it is customary 
to speak about gauge dependence of the Wannier functions due to the /c-dependence of 
the U{n) transformations. In the single band case the arbitrariness reduces to the 
abelian f/(l) group of phase transformations coming from the freedom in choosing the 
Bloch basis at each k. 

A general approach to obtain MLWFs has been proposed in a seminal paper by 
Marzari and Vanderbilt ^21j , where MLWFs are obtained by minimizing the generalized 
Wannier spread VL = [{x'^)n — {x)n] by means of a suitable gauge transformation of 
the Bloch eigenf unctions. In the case of a single band the method returns the Kohn 
result [21] ■ Marzari and Vanderbilt have shown that the spread can be written as the 
sum of two positive terms, Q = Qj + Q, the first being gauge invariant and therefore 
fixing the minimal spread. The gauge dependent term Q can be further split into the 




(20) 




(21) 



(22) 
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diagonal and off-diagonal components, Q = Qo + ^od, that in the one- dimensional case 
read 

= ^^\{Wnj\x\Wno) (23) 

n j^O 

^OD = ^^^{Wmj\x\Wno) ■ (24) 
mj^n j 

Both Qj:, and Qqd can be written in terms of the generalized Berry vector potentials 
Anm{k), defined as [32| 133] 

Anm{k) = l—{Unk\dk\Umk) (25) 

with the matrix A{k) being hermitian, A'^lk) = A{k) (it follows from dk{unk\um,k) = 0). 
We also recall that the integral over the first Brillouin zone of Ann{k) gives the one band 
Zak- Berry phases 7^ 



7^ 



n = i—r I {Unk\dk\Unk) = / Ann{k) = -^(^nn)^ (26) 

d Jb Jb d 

which are proportional to the offset of the Wannier function centers, (x)„o = 
{{x)nj — Rj) = {d/27i)'jn [311 El], yielding {x)no = {Ann)B (this relation is preserved 
under a generic unitary gauge transformation, as in (1221) ). It is also worth to remember 
that the single Wannier centers are invariant only under single band U{1) gauge 
transformations, whereas for general U (n) transformations only their sum is conserved 

Then, the expressions for and Qqd can be written as 

nn = J2((^nn{k) - {Ann}Bf}B = J^^^n (27) 
n n 

^OD = (l^nmHa , (28) 

and in general can be reduced by means of a functional minimization in k space, as 
discussed in [2ll|22]. 

Here we use a different approach, specifically suited for constructing the set of 
MLWFs for the double-well case. Namely, we show that the minimization problem can 
be reformulated by identifying a specific gauge transformation for a composite band 
in one dimension, expressed in terms of a set of ordinary differential equations with 
periodic boundary conditions. We recall that in one dimension Cl can be made strictly 
vanishing, and this corresponds to find a gauge (also called parallel transport gauge) in 
which the matrix Anm{k) is diagonal, with the diagonal elements being constant and 
equal to their mean values. The latter are related to the eigenvalues of the matrix 
generalizing the Berry phase to the non abelian case [21] . 

3.1. Gauge transformations and differential equations 

The diagonal and off-diagonal spreads Qd and Qqd can be minimized either 
simultaneously or independently. For the following discussion, it is useful to distinguish 
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between two kind of gauge transformations. 

I. Diagonal U{n) transformations which correspond to a set of single band gauge 
transformations of the form 

\Unk) -> \Unk)=e''^"^''^\Unk) (29) 

with (j)n{k) being a real continuous (differentiable) function of k, such that (f)n{k + 
2^5) = <Pn{k) + 27r£ {[ integer) in order to have periodic and single valued Bloch 
eigenstates. Then, we can set £ = without loss of generality. As discussed in 
[2T] . this transformation may be used to minimize each term of VLd, as it affects 
the Wannier spread il„ = — (x),^ while preserving their centers (x)„ (modulo 

a lattice vector). In particular, VLo can be set exactly to zero [21]. In fact, we have 
Ann{k) Ann{k) - dk(f>n{k), and therefore 

Qon ^Dn = {{Ann - dk4>n - {Ann) Bf) B (30) 

that vanishes by imposing 

9fc0„ = Ann - {Ann)B ■ (31) 



This equation can be readily solved numerically, as discussed in [Appendix A| and 



Appendix B In addition, it is straightforward to verify that under the transformation 
(|29|) . Ai2{k) changes only by a phase factor, and therefore VLod remains unchanged. 
II. A full gauge transformation in the composite band of the form 

\Unk) \Unk) = ^Unm{k)\Unik) (32) 

m 

where, as already said, U{k) G U{N) (for a A^-composite band), constrained to the 
following periodic condition 

Unm {k + 2kB) = Unm{k) . (33) 

Under such a transformation the generalized Berry potentials transform as 

Anm ^ Anm i / dx Un^kl^mk 

a J 

= I KidkUn^i + KlUmvAu, (34) 
/ 1,1' 

In general, U{N) can be written as a semidirect product SU{N) x f/(l), with 
U{1) being subgroup of U{N) consisting of matrices of the form diag(l, 1, . . . , e''^). As 
anticipated, here we restrict to a 2 x 2 case, and therefore we have 

with + |z3p = 1, r{k) = e^^^^\ Moreover, by using the following parametrization 
for a matrix S G SU{2), S = e^^^ ' with n = (cos sin6', siny^ sin 6', cos^) and ai 
being the Pauli matrices, we can write 

cos I + « sin I cos 6* ze^(^ ~ sin 6* sin | | ^^^^ 
\ ie^*^ sin 9 sin | e*^ (cos | — i sin | cos 9) J 
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with X = x{k), = V^(^), a = o:{k) and 6 = 6{k). 
Then, since Q transforms as follows 

^Dn ^Dn = ((^nn(A;) - (inn)^)^)^ (37) 

^ = 2{\Au\^)s (38) 
in order to get 17 = one has to impose (see ([21])) 

Annik) = I J2 KldkUnl + K.Unl' A, = (Ann) B (^ = 1, 2) (39) 

I 1,1' 

Au{k) = I J2 UIAU21 + U*iU2i'Au' = 0. (40) 
I 1,1' 

Notice that the former, (139|) . is an integro-differential equation where the right-end side 
corresponds to the center of the Wannier functions, (A„„)g = {wnj\x\wnj) , that are not 
known a priori (only their sum is conserved in the parallel transport gauge). 

Therefore, in the following we will consider a specific transformation that makes 
Qqd vanish without specific requirements on the transformation of Qjj, given by the 
solution of (HOl) . Then, by using (136|) . the latter can be transformed in a system of four 
differential equations for a, 6, ip, x, whose normal form is 

dk(y cos 26 , .r, .j. x 

~ cos r] + sm r/) 

- cotg^cotg6' (y4f2 sin rj - cos rj) + cos 9{Aii - A22) (41) 

„ „ cos sin a; . ^ p . , 

= . 21 ,n\ (^12 cos T] + A-^2 sm ri) 
sm (a/2) 

+ . 2( ,nM i2 smr/ - cosr^) - cotg- sm^(An - A22) (42) 
sm (a/2) 2 

9fcX = 0, dkV = Q (43) 

where we have defined rj = ip — x, with dkri = 0. The solution of is x = Xoy 
if = ifo. Then, it is evident that only two equations are left, namely fHT]) and f l42|) . with 
rj = (po — Xo playing the role of a parameter, and with one of the angles among (po and 
Xo being arbitrary. We then can choose Xo = without loss of generality. In addition, 
in order to conform to ( 133|) . the angles a and 9 must satisfy the following periodicity 
conditions 

a{-kB) = a{kB) + Mtt (44) 
Oi-ks) = e{kB) + 2f TT (45) 

with £, ^' integers which can be taken £ = = without loss of generality. This set 
of equations, fllT]) . fH2]) . (jH]), and can be solved as discussed in Appendix B The 
resulting transformation is an SU{2) matrix S{k) of the form fl36|) with x = a-^id ^9 
constant. 

Summarizing, the procedure to make Q vanish can be divided in two steps: (a) a. 
gauge transformation of type II to make Qqd vanish (without specific requirements 
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Figure 3. Plot of the density of the two lowest single band (blue lines) and generalized 
(red lines) MLWFs, in log (a,b) and linear scale (c,d). The dotted line in (c,d) represent 
the potential, while the horizontal orange stripes are the first three Bloch bands (on 
the same scale of the potential). Here Vi — IQEj^, V2 = 20i?/f, 9o = Tr/2. 



on the transformation of the diagonal elements Ann{k)); (h) a set of two gauge 
transformation of type I (one for each band) that makes vanish without affecting 
VLqd- Therefore, the full transformation can be decomposed in the following way 

U^m{k) = e''^-^^)s^m{k) (46) 



It is straightforward to verify that fH6l) can be cast again in the form fl35|) by properly 
redefining the parameters 2:1, 2:3, r. 

3.2. Examples of MLWFs 

Here we present some example of generalized MLWFs obtained with the transformation 
discussed above, and we compare them with the single hand MLWFs. The latter can 
be obtained by using just the gauge transformation of type I, and correspond to the 
exponentially decaying Wannier functions discussed by Kohn p!3| [2T] . Both gauge 
transformations are solved by using the representation of Bloch functions in fc-space, by 



means of the numerical methods discussed in Appendix A and Appendix B 

In figure Owe show the case of a symmetric double well, 6q = 71 /2. In this case the 
single band MLWFs have the same symmetry of the potential [13], therefore occupying 
both wells in the unit cell. Instead, each of the generalized MLWFs nicely localizes in 
one of the sub-wells, as wanted 



§ We remark that in this symmetric case one could construct a set of localized functions resembling 
the generalized MLWFs by considering symmetric and antisymmetric combinations of the single band 
MLWFs. These combinations (not shown) nicely reproduce the bulk density of the generalized MLWFs 
when each unit cell can be regarded as a single double- well (that is, large barriers at the cell borders). 
However, this approximation does not work in the region of the tails (those of the generalized MLWFs 
decaying much faster), therefore not being reliable for computing the tunnelling coefficients. It is also 
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Figure 5. Same of figure [3] for Oq = 0,Vi = lEn. 



Then, in figures IH |5] we show two cases of an asymmetric double well, here for 
Oq = 0. These figures are useful for illustrating the role of the band gaps. The case in 
figure m corresponds to a band gap between the first two bands that is larger than the 
one between the second and third band; in this case both the single band MLWFs are 
already localized within a single well, the one in panels (b,d) having a small residual 
amplitude around the neighbouring wells. The generalized MLWFs do not differ much 
from the former in linear scale, the effect of the gauge mixing transformation being a 
reduction of the lateral lobes of the single band MLWF in (b,d), but the price to pay is 
a substantial increase of the width of the other one in log scale, see panel (a). Actually, 
in this case composite band approach is not fully justified, due to the large band gap 

possible to prove analytically that this approximation implies Jab~ — Jab+, that is manifestly in 
contradiction with the definition of the tunneling coefficients given in figure [H Nevertheless, such wave 
functions may be useful for discussing the boundary conditions for (|¥T]) - (P5)) . see [Appendix B 
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Figure 6. Density plot of tlie ratio between tlie second and first band gaps, 
R = £gi2/£g23 a function of 9 and V2, for Vi = bEn (a), and Vi = lOEji (b). 
The black line corresponds to i? = 1. The colour scale is saturated at i? = 1. 

between the first and second band, and one should consider the structure of the upper 
bands. This however goes beyond the scope of this work. 

Instead, figure [5] shows a case still for 9o = but where an almost degeneracy 
between the two minima (that corresponds to a quasi degeneracy of the lowest two 
bands) is restored thanks to a lower value of Vi (here Vi = 1 instead of 10 as in the 
former figure). Here the advantage of using the generalized MLWFs is clearly visible. 

4. Tight binding regimes and tunnelling coefficients 

In this section we will discuss the features of the composite band approach, by comparing 
the predictions of the full and nearest-neighbour versions of the model, and discussing 
the behaviour of the tunnelling coefficients and the other tight-binding parameters. 

We recall that the potential is characterized by three parameters, Vi/Er, V2/ER 
and 6*0. The latter can be restricted to the range [0,7r/2] without loss of generality, see 
figured! As for the potential amplitudes, we chose 10 < V2/ Eji < 40 in order to fulfil the 
tight binding regime, in a feasible range for current experiments with ultra cold atoms 
[1]. In addition, for having both wells deep enough in order to have the corresponding 
MLWFs mostly localized within each well, Vi should not exceed V2/2. 

As we have anticipated in the previous section, in principle the composite band 
approach is justified in a situation of quasi degeneracy between the first two Bloch 
bands, that is when their separation Sg^^ is much smaller that the band-gap £^23 between 
the second and third band. Therefore, it is convenient to define the ratio R = £312/^923) 
whose behaviour is shown as a density plot in figure [6] as a function of 6q and V2, for 
different values of Vi. As a matter of fact, we find that the composite band approach 
provides a basis of functions well localized in each of the two sub-wells and correctly 
reproduces the lowest two Bloch energy bands even up to i? ~ 1, in the parameter range 
discussed above. 

In order to characterize the level of fidelity in reproducing the exact single particle 
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Figure 7. Plot of the quantity fc„ (see text) for Vi = 5, as a function of for 
V2 = 20 (a), and as a function of V2 ioT 9q — tt/2 (b). Red symbols: full tight-binding 
model; Blue symbols: nearest-neighbour approximation. Empty squares: first band; 
Solid squares: second band. These figures refer to a horizontal and vertical cut in the 
left panel of figure ID 



Bloch spectrum (that can be readily computed as discussed in Appendix A ), we define 
the following quantity 




^^n^-r-xhz dk[eM-et{k)f (47) 



that represents the ratio of the quadratic spread between the exact Bloch spectrum 
Snik) and that of the tight binding hamiltonians ([8]) and ([9]), to the bandwidth 
Ae^ = (e™"^ — e™"). We also notice that the formula f fTTl) for e^ik) provides the 
same numerical result of the expression (1C.6P for the energy bands in terms of the 
gauge transformations, and that these formulas have a better accuracy in reproducing 
the exact Bloch spectrum than the single band expression (ITTl) . especially in the region 
close to the symmetric case = 7r/2. 

The quantity Ssn is shown in figure[7]for Vi = 5, as a function of ^0 = 20) and V2 
(^0 = ^/2), in the left and right panels respectively. These figures refer to a horizontal 
and a vertical cut in the left panel of figure [6l and confirm that in the regime i? < 1 
the full tight binding model reproduces the exact energies with great accuracy. As for 
the nearest neighbour approximation commonly used in the literature, in general this is 
not the it works reasonably only for i? < 0.1, that is for 9q ~ 7r/2 and large V2. 

Therefore, attention must be payed to the regimes where it is allowed to neglect some 
of the next-to-nearest tunnelling terms. 



Tunnelling coefficients 

Let us now consider the explicit expressions for the onsite energies and the tunnelling 
coefficients in ([8]). They can be expressed in terms of the gauge transformations discussed 
in the previous section as (see and (jlSD) 

d r ^ 

=7r dky2\S,Uk)\'em{k) (48) 



Maximally localized Wannier functions for ID double-well periodic potentials 15 




0.1 0.2 0.3 0.4 0.5 10 20 y 30 40 



Figure 8. Plot of the tunnelling coefficients (in modulus, rescaled to Tab) for Vi = 5, 
as a function of 6*0 for V2 = 20 (a), and as a function of V2 for — tt/2 (b). These 
figures refer to a horizontal and vertical cut in the left panel of figure [5] 




Figure 9. Plot of the modulus of Tab (xlOO) and of the onsite energy difference 
Eb — Ea, for the same parameters of the previous figure. 



J. 

Tab 
Jab± 

with u = A, B = 1,2 and A(f){k) = (j)2{k) — 4>i{k). We recall that all these terms are 
relevant for the full tight-binding model, while the nearest-neighbour approximation 
commonly used in the literature corresponds to retain just the contribution of E^, Tab 
and Jab_- We also notice that all the above parameters are gauge dependent; and 
Jy depend just on the gauge mixing transformation; Tab and Jab± on the whole gauge 
transformation. 

It is worth mentioning that in the parallel transport gauge only the modulus of 
the tunnelling coefficients is univocally defined. However, as it is briefiy discussed in 
Appendix B, there is the freedom to choose them real. 

The behaviour of the tunnelling coefficients (in modulus, rescaled to Tab) is shown 
in figure [S] as a function of 9q and V2 (for V2 = 20 and 9q = 7i/2, respectively; here 



^ [ dke-''''y^\SUk)?em{k) (49) 
— / dk e'^'^'^'^ V SUk)S2m{k)em{k) (50) 
I g.(A0(/c)Tfcd) S*Uk)S2^^{k)e^{k) (51) 



d_ 
2^ 
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9/71 

Figure 10. Plot of the rescaled interaction term I^, for Vi = 5, as a function of 6*0 
for V2 — 20 (a), and as a function of V2 for Oq — tt/2 (h). Empty squares: first band; 
Solid squares: second band. In the right panel the two lines are superimposed owing 
to the double well symmetry. 



Vi = 5), whereas the absolute variation of Tab and of the onsite energy difference 
Eb — Ea is shown in figure O The former figure reveals that the dependence on ^0 is 
rather weak, the only notable effect being that at ^0 = (where all the maxima are 
degenerate) Tab = Jab+, whereas at 9q = Tr/2 we have Ja = Jb (and Ea = Eb, as in 
this case the minima are degenerate). Instead, the tunnelling coefficients Ja, Jb and 
Jab+ are substantially suppressed with respect to Tab and Jab- as V2 is increased, and 
consequently the results of the nearest-neighbour approximation improve. Also, figure 
[7^ shows that the bands energies obtained with the nearest-neighbour approximation 
better reproduce the exact results by approaching 6*0 = 7r/2 in spite of the very weak 
variation of all the tunnelling coefficients shown in figure 13 Actually, a detailed analysis 
reveals that these results are due to the reduction of the band gap Eg^^ and to the increase 
of the lowest bands widths A^i 2 when going from 6'o = towards 9q = 7^/2. Both from 
(|T7|) and (1C.6I) it is possible to realize that the relative weight of the terms containing 
J A and Jb becomes consequently less relevant. 

Finally, in figure fTOl we consider the behaviour of the rescaled interaction term (see 
ffTT]) ). Ii, = U^/g = j dx \wj^{x)\'^, again as a function of 6q and V2, for the parameters 
of the preceding figures. This term is also called inverse participation ratio (IPR), and 
represent the number of occupied lattice sites in the case of a model defined on a discrete 
lattice (see e.g. [71I3B]). In the present continuous case, values of /j, close or greater than 

I imply that the MLWFs are localized on a distance smaller than the lattice spacing. 

5. Conclusions 

We have presented a method for constructing a set of maximally localized Wannier 
functions (MLWFs) for a one- dimensional periodic potential with a double-well per 
unit cell. Starting from the minimal spread requirement by Marzari et al. |21j . 

II In the regime of parameter considered here, the next-to- leading interaction terms - see pop - are 
suppressed by about 2 orders of magnitude. 
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we have designed a two-step gauge transformation specific for a composite two band 
system, consisting a set of ordinary differential equations with periodic boundary 
conditions. In the tight binding regime, this method provides a proper set of MLWFs 
when both wells are sufficiently deep, provided that the band gap between the first 
two bands is larger than the one between the second and third band. By using 
these MLWFs one can map the continuous hamiltonian onto a tight binding model, 
whose coefficients are expressed in terms of the above gauge transformations, providing 
a direct way to relate the tight binding coefficients to the experimental parameters 
characterizing the continuous potential. We have also shown that in general, in the 
range of typical experimental parameters, it is necessary to retain all possible tunnelling 
couplings between the different sub- wells of neighbouring cells [full model), whereas 
the usual nearest neighbouring approximation used in the literature is justified only for 
V2/E11, V2/V1 3> 1 and 9o ~ 7i/2. The method can be applied also in higher dimensions, 
in case of separable lattice potentials. 
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Appendix A. Numerics: discretization in momentum space 

In order to address the problem from the numerical point of view, we use the following 
truncated Fourier expansion of the Bloch eigenf unctions 

Un,ix) = ^Y1 (A.l) 



2^ e=-N 

with J2e = 1- The eigenvalue equation H\unk) = £n{k)\unk) is then mapped onto 
a matrix eigenvalue equation for the coefficients Cne{k), and can be solved by means of 
standard linear algebra routines, like LAPACK. Then, the Berry connections can be 
written as 

and can be computed on a /c-grid with spacing Ak 0, by means of symmetric two-point 
derivative 

dcmi{k) Cme{k + Ak) - Cme{k - Ak) 

" 2Ak (^-'^ 

by using the periodicity condition Cni{k + 2kB) = c„£+i(fc), that follows from the fact 
that the ipnk{x) are periodic in quasimomentum with period 2k b, or equivalently that 

Unk+2kB{x) = e-'^''^''Unk{x) (A.4) 



% Typically have used from 500 to 10000 grid points. 
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Since the diagonalization at each point in fc-space is uncorrelated from that at 
neighbouring ks, in general the diagonahzation routines return a set of coefficients Cme{k) 
with an arbitrary phase factor e*'''"'" that cannot be controlled a priori, and that affects 
their differentiability in (IA.3p . However, a set of smooth (differentiable) Cmiik) can be 
obtained as follows. We notice that in a; = the periodicity condition f lA.4|) reads 
^n/c+2fcs(0) = Uk{0). Therefore, given a set of Bloch eigenfunctions m„a;(x) generated 
numerically (affected by the same spurious phase factors e*'^"'= of the c„^(/c), see flA.ip ) 
it is possible to define a new set of smooth Bloch eigenfunctions Unk{x) by fixing all the 
phases to zero at x = (i.e. by imposing them to be real in x = 0), 

Unk{x) = Unk{x) exp(-i arg(M„fc(0))) (A. 5) 

therefore fixing the spurious phases. Formally, this corresponds to a diagonal gauge 
transformation generated by (pnik) = — arg(M„fc(0)) = — arg (^^ c„f (A;)). Since 0„(/c) 
does not depend on i, the Cni{k) transform in the same way 

Cni{k) = Cni{k) exp ^-i arg c„^(A:)^ j (A.6) 

and this can be used as a prescription for ensuring the smoothness of the Cne{k) and Unk 
as a function of k. 



Appendix B. Numerics: gauge transformations 

Diagonal transformation. The diagonal gauge transformation of type I is obtained 
by solving the ordinary differential equation for (pnik) in fl3T]) by using a 4th-order 
Runge-Kutta (RK4) integrator with initial condition (f)n{—kB) = 0. A single run 
returns a (single band) diagonal spread flDn = ^ > 0, with e ^ 1; the procedure can 
be iterated in order to achieve the desired accuracy e^m, ^Dn < ^min (typical value: 
e^,„ = 10-18 ^10-25). 

Gauge mixing transformation: method 1. The solution of fHTl) and fH2|) with 
periodic boundary conditions dS]), (jlSD (here we fix ^ = ^' = 0) can be achieved 
by using two nested shooting algorithms, adapting the discussion in [35j to the case 
considered here. As a first step ("inner loop"), we fix a trial value for 77, and we chose 
an arbitrary initial value a(— /c^) = olq (avoiding the singular value ao 7^ ^r). Then we 
impose the periodic boundary conditions a(— fc^) = a{kB) with a shooting method that 
uses the initial value 9{—kB) as free parameter, that is varied in order to match the 
required boundary conditions on a [35]. The differential equations are integrated with 
the same RK4 integrator mentioned above. The second step ("outer loop") consists 
in a similar shooting algorithm, that uses r/ as a free parameter trying to match the 
periodic boundary conditions on 6, O^—ks) = O^ks), without affecting those on a (that 
are constrained by the "inner loop"). With this method it is possible to reduce Qqd 
by several orders of magnitudes, down to ^ 10"^^, (the actual value depends on the 
number of grid points and on the region of parameters). We also find that the method 
is more efficient if it is preceded by the diagonal transformation discussed above. 
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Gauge mixing transformation: method 2. In order to solve (1411) . fH2l) . we have 
also employed another method that originates from the observation that in the case of 
two degenerate symmetric minima in the elementary cell (^o = 7r/2), simple symmetric 
and antisymmetric combinations of Bloch eigenstates give rise to two Wannier functions 
located each at one of the two degenerate minima (though they are not the maximally 
localized ones). In this situation, we have verified that initial conditions for the matrix 
Snmi^kB) - corresponding to symmetric and antisymmetric combinations of Bloch 
eigenstates - allow to obtain the periodic solutions minimizing Qqd which coincide 
both integrating from —ks to ks and vice- versa. By smoothly changing 6q, we start 
from the same initial conditions and then we iterate the integration of fHTj) . fH2|) by 
taking as initial angles, at each step, the mean values of the final values of the preceding 
integration (at the two boundaries of the BZ). The procedure rapidly converges with 
QoD reduced by several order of magnitude (this depends on the number of grid points 
and on the region of parameters, as for the previous method). In particular, for 6q = 7r/2 
(symmetric double well) we have always found solutions with ipo = n/2 and the following 
boundary conditions for a and 6 

a{±kB) = e{±kB) = n/2. (B.l) 

By smoothly changing ^o, the corresponding boundary conditions are smoothly 
connected with this set (in the range of Vi and V2 that we have explored). Notably, 
this choice yields real tunnelling coefficients (in principle this is not guaranteed, see 
Appendix 0.2]). 



Appendix C. Uniqueness of the physical results 
Appendix C.l. Energy bands in the truncated expansion 

Let us start by considering the following exact expression for the Bloch spectrum 

e^{k) = J2e"'''Ue) (C.l) 

e 

with the coefficients T„(£) being the expectation values of on single band Wannier 
functions [H] 

T„(£) = («;„,|i^okn,+£) = ^ [dk e^{k)e-'^^'. (C.2) 

2vr Jb 

This expression is manifestly gauge invariant. Instead, the analogous expression in terms 
of the generalized MLWFs depends on the gauge. In fact, with some manipulations we 
can write (here we are implicitly considering just the first two bands, n = 1,2) 

£„(A;) = ^e*^^'^P„(A;,£) (C.3) 

with 

Pn{k,l) = Y,Un"n{k)U:,^{k)fn"n'{i) (C.4) 
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and 

Tn"n'{^) = {Wn''j\Ho\Wn'j+e). (C.5) 

By truncating (1C.3P to nearest neighbouring cells {i = 0, ±1) and using the definitions 
of the onsite energies and of the tunnelling coefficients in ©-([T]), we have 

Bnik) = |f/inP^A + |f/2n|'^B " 2\Uin\^Re [c'^'' J a) - 2|f/2„pRe (e^'^Vs) 

- 2Re [f/i„f/*„ [Tab + e^'V^B^ + e-^^^V^fiJ] (C.6) 
with Unm{k) = e'*^-^^^ Snm{k) and S e SU{2) (see Sect. [3]). Though each energy 6n{k) 
in (IC.6P depends on the gauge (through the the matrix Unm{k)), it is not difficult 
to demonstrate that their sum, J2n=i 2^n{k), is gauge invariant like the analogous 
expression obtained form (IC.ip . In addition, once we have chosen the parallel transport 
gauge, corresponding to the generalized MLWFs with minimal spread {Q = 0), the 
expression f lC.6p is univocally defined despite the presence of residual phase arbitrariness 
(see next paragraph). 



Appendix C.2. Residual phase ambiguities 

Let us now show that the relevant quantities are univocally defined in the parallel 
transport gauge, despite possible ambiguities in the procedure for obtaining the matrix 
Unm{k). First of all, we remind that the whole procedure starts from a particular gauge. 



as discussed in Appendix A Choosing a different initial gauge amounts to perform a 
diagonal transformation on Bloch eigenstates 

\i'mk) ^ e'^-^^^i^mk) (C.7) 

being (pmik) an arbitrary periodic real function (see Sect. [3]). Consequently, the full 
gauge transformation matrix leading to the parallel transport gauge changes as 

Unm{k) ^ f/„„(A;)e-^*™('=). (C.8) 

Then, it is easy to see that the onsite energies and tunnelling coefficients in ([S])- 
([7]) are invariant under this transformation, as well as the energies in fIC.Gp . whereas 
the real and imaginary parts of the MLWFs change, conserving their modulus. The 
procedure to make VL vanishing does not imply any other fc-dependent ambiguity on the 
matrix Unm{k)- In addition, in the parallel transport gauge we still have the freedom 
to perform two independent U{1) gauge transformations for each band with angles 
which necessarily obey dXn/dk = 0. These terms appear as multiplicative factors on the 
rows of the matrix Unm- Finally, it turns out that - given a set of parameters Vi, V2 and 
^0 - there is not a unique set of boundary conditions for the angles (fo, « and 6 satisfying 
( HTj) . f H2]) (and thus not an unique transformation Snmi^ks) leading to ^od = 0). For 



^0 = 7r/2 these boundary values can be fixed as discussed in Appendix B Different 
choices may also be possible. However, it is not difficult to prove that none of these 
ambiguities affect either the energies in f lC.6p . or the modulus of the MLWFs and of the 
tunnelling coefficients. Instead, their real and imaginary parts are not univocally defined 
and can be chosen somewhat arbitrarily. In particular, real tunnelling coefficients can 



be obtained by using the set of boundary conditions given at the end of Appendix B 
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